// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#include <stdlib.h>
#include <string.h>
#include <map>
#include <string>
#include <sstream>
#include <cassert>
#include <iomanip>
#include "GModel.h"
#include "OS.h"
#include "GmshDefines.h"
#include "MPoint.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "MElementCut.h"
#include "StringUtils.h"
#include "GmshMessage.h"
#include "Context.h"
#include "OS.h"

// periodic nodes and entities backported from MSH3 format
extern void writeMSHPeriodicNodes(FILE *fp, std::vector<GEntity *> &entities,
                                  bool renumber);
extern void readMSHPeriodicNodes(FILE *fp, GModel *gm);
extern void writeMSHEntities(FILE *fp, GModel *gm);

static bool getMeshVertices(int num, int *indices,
                            std::map<int, MVertex *> &map,
                            std::vector<MVertex *> &vertices)
{
    for(int i = 0; i < num; i++) {
        if(!map.count(indices[i])) {
            Msg::Error("Wrong vertex index %d", indices[i]);
            return false;
        }
        else
            vertices.push_back(map[indices[i]]);
    }
    return true;
}

static bool getMeshVertices(int num, int *indices, std::vector<MVertex *> &vec,
                            std::vector<MVertex *> &vertices, int minVertex = 0)
{
    for(int i = 0; i < num; i++) {
        if(indices[i] < minVertex ||
                indices[i] > (int)(vec.size() - 1 + minVertex)) {
            Msg::Error("Wrong vertex index %d", indices[i]);
            return false;
        }
        else
            vertices.push_back(vec[indices[i]]);
    }
    return true;
}

static MElement *
createElementMSH2(GModel *m, int num, int typeMSH, int physical, int reg,
                  unsigned int part, std::vector<MVertex *> &v,
                  std::map<int, std::vector<MElement *> > elements[10],
std::map<int, std::map<int, std::string> > physicals[4],
bool owner = false, MElement *parent = 0, MElement *d1 = 0,
MElement *d2 = 0)
{
    if(CTX::instance()->mesh.switchElementTags) {
        int tmp = reg;
        reg = physical;
        physical = tmp;
    }

    MElementFactory factory;
    MElement *e = factory.create(typeMSH, v, num, part, owner, 0, parent, d1, d2);

    if(!e) {
        Msg::Error("Unknown type of element %d", typeMSH);
        return nullptr;
    }

    int dim = e->getDim();
    if(physical &&
            (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
        physicals[dim][reg][physical] = "unnamed";

    if(part > m->getNumPartitions()) m->setNumPartitions(part);
    return e;
}

int GModel::_readMSH2(const std::string &name)
{
    FILE *fp = Fopen(name.c_str(), "rb");
    if(!fp) {
        Msg::Error("Unable to open file '%s'", name.c_str());
        return 0;
    }

    char str[256] = "XXX";
    double version = 1.0;
    bool binary = false, swap = false, postpro = false;
    std::map<int, std::vector<MElement *> > elements[10];
    std::map<int, std::map<int, std::string> > physicals[4];
    std::map<int, MVertex *> vertexMap;
    std::vector<MVertex *> vertexVector;
    int minVertex = 0;

    while(1) {
        while(str[0] != '$') {
            if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
        }

        if(feof(fp)) break;

        if(!strncmp(&str[1], "MeshFormat", 10)) {
            if(!fgets(str, sizeof(str), fp)) {
                fclose(fp);
                return 0;
            }
            int format, size;
            if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) {
                fclose(fp);
                return 0;
            }
            if(format) {
                binary = true;
                Msg::Debug("Mesh is in binary format");
                int one;
                if(fread(&one, sizeof(int), 1, fp) != 1) {
                    fclose(fp);
                    return 0;
                }
                if(one != 1) {
                    swap = true;
                    Msg::Debug("Swapping bytes from binary file");
                }
            }
        }
        else if(!strncmp(&str[1], "PhysicalNames", 13)) {
            if(!fgets(str, sizeof(str), fp)) {
                fclose(fp);
                return 0;
            }
            int numNames;
            if(sscanf(str, "%d", &numNames) != 1) {
                fclose(fp);
                return 0;
            }
            for(int i = 0; i < numNames; i++) {
                int dim = -1, num;
                if(version > 2.0) {
                    if(fscanf(fp, "%d", &dim) != 1) {
                        fclose(fp);
                        return 0;
                    }
                }
                if(fscanf(fp, "%d", &num) != 1) {
                    fclose(fp);
                    return 0;
                }
                if(!fgets(str, sizeof(str), fp)) {
                    fclose(fp);
                    return 0;
                }
                std::string name = ExtractDoubleQuotedString(str, 256);
                if(name.size()) setPhysicalName(name, dim, num);
            }
        }
        else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5) ||
                !strncmp(&str[1], "ParametricNodes", 15)) {
            const bool parametric = !strncmp(&str[1], "ParametricNodes", 15);
            if(!fgets(str, sizeof(str), fp)) {
                fclose(fp);
                return 0;
            }
            int numVertices = -1;
            if(sscanf(str, "%d", &numVertices) != 1) {
                fclose(fp);
                return 0;
            }
            Msg::Info("%d vertices", numVertices);
            Msg::ResetProgressMeter();
            vertexVector.clear();
            vertexMap.clear();
            minVertex = numVertices + 1;
            int maxVertex = -1;
            for(int i = 0; i < numVertices; i++) {
                int num;
                double xyz[3], uv[2];
                MVertex *newVertex = 0;
                if(!parametric) {
                    if(!binary) {
                        if(fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) !=
                                4) {
                            fclose(fp);
                            return 0;
                        }
                    }
                    else {
                        if(fread(&num, sizeof(int), 1, fp) != 1) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)&num, sizeof(int), 1);
                        if(fread(xyz, sizeof(double), 3, fp) != 3) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
                    }
                    newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
                }
                else {
                    int iClasDim, iClasTag;
                    if(!binary) {
                        if(fscanf(fp, "%d %lf %lf %lf %d %d", &num, &xyz[0], &xyz[1],
                                  &xyz[2], &iClasDim, &iClasTag) != 6) {
                            fclose(fp);
                            return 0;
                        }
                    }
                    else {
                        if(fread(&num, sizeof(int), 1, fp) != 1) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)&num, sizeof(int), 1);
                        if(fread(xyz, sizeof(double), 3, fp) != 3) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
                        if(fread(&iClasDim, sizeof(int), 1, fp) != 1) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)&iClasDim, sizeof(int), 1);
                        if(fread(&iClasTag, sizeof(int), 1, fp) != 1) {
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)&iClasTag, sizeof(int), 1);
                    }
                    if(iClasDim == 0) {
                        GVertex *gv = getVertexByTag(iClasTag);
                        if(gv) gv->mesh_vertices.clear();
                        newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
                    }
                    else if(iClasDim == 1) {
                        GEdge *ge = getEdgeByTag(iClasTag);
                        if(!binary) {
                            if(fscanf(fp, "%lf", &uv[0]) != 1) {
                                fclose(fp);
                                return 0;
                            }
                        }
                        else {
                            if(fread(uv, sizeof(double), 1, fp) != 1) {
                                fclose(fp);
                                return 0;
                            }
                            if(swap) SwapBytes((char *)uv, sizeof(double), 1);
                        }
                        newVertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, uv[0], num);
                    }
                    else if(iClasDim == 2) {
                        GFace *gf = getFaceByTag(iClasTag);
                        if(!binary) {
                            if(fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) {
                                fclose(fp);
                                return 0;
                            }
                        }
                        else {
                            if(fread(uv, sizeof(double), 2, fp) != 2) {
                                fclose(fp);
                                return 0;
                            }
                            if(swap) SwapBytes((char *)uv, sizeof(double), 2);
                        }
                        newVertex =
                                new MFaceVertex(xyz[0], xyz[1], xyz[2], gf, uv[0], uv[1], num);
                    }
                    else if(iClasDim == 3) {
                        GRegion *gr = getRegionByTag(iClasTag);
                        newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
                    }
                }
                minVertex = std::min(minVertex, num);
                maxVertex = std::max(maxVertex, num);
                if(vertexMap.count(num))
                    Msg::Warning("Skipping duplicate vertex %d", num);
                vertexMap[num] = newVertex;
                if(numVertices > 100000)
                    Msg::ProgressMeter(i + 1, numVertices, true, "Reading nodes");
            }
            // If the vertex numbering is dense, transfer the map into a
            // vector to speed up element creation
            if((int)vertexMap.size() == numVertices &&
                    ((minVertex == 1 && maxVertex == numVertices) ||
                     (minVertex == 0 && maxVertex == numVertices - 1))) {
                Msg::Debug("Vertex numbering is dense");
                vertexVector.resize(vertexMap.size() + 1);
                if(minVertex == 1)
                    vertexVector[0] = 0;
                else
                    vertexVector[numVertices] = 0;
                std::map<int, MVertex *>::const_iterator it = vertexMap.begin();
                for(; it != vertexMap.end(); ++it) vertexVector[it->first] = it->second;
                vertexMap.clear();
            }
        }
        else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) {
            if(!fgets(str, sizeof(str), fp)) {
                fclose(fp);
                return 0;
            }
            int numElements;

            std::map<int, MElement *> elems;
            std::map<int, MElement *> parents;
            std::map<int, MElement *> domains;
            std::map<int, int> elemreg;
            std::map<int, int> elemphy;

            std::set<MElement *> parentsOwned;
            sscanf(str, "%d", &numElements);
            Msg::Info("%d elements", numElements);
            Msg::ResetProgressMeter();
            if(!binary) {
                for(int i = 0; i < numElements; i++) {
                    int num, type, physical = 0, elementary = 0, partition = 0,
                            parent = 0;
                    int dom1 = 0, dom2 = 0, numVertices;
                    std::vector<short> ghosts;
                    if(version <= 1.0) {
                        if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary,
                                  &numVertices) != 5) {
                            fclose(fp);
                            return 0;
                        }
                        if(numVertices != (int)MElement::getInfoMSH(type)) {
                            fclose(fp);
                            return 0;
                        }
                    }
                    else {
                        int numTags;
                        if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) {
                            fclose(fp);
                            return 0;
                        }
                        int numPartitions = 0;
                        for(int j = 0; j < numTags; j++) {
                            int tag;
                            if(fscanf(fp, "%d", &tag) != 1) {
                                fclose(fp);
                                return 0;
                            }
                            if(j == 0)
                                physical = tag;
                            else if(j == 1)
                                elementary = tag;
                            else if(version < 2.2 && j == 2)
                                partition = tag;
                            else if(version >= 2.2 && j == 2 && numTags > 3)
                                numPartitions = tag;
                            else if(version >= 2.2 && j == 3)
                                partition = tag;
                            else if(j >= 4 && j < 4 + numPartitions - 1)
                                ghosts.push_back(-tag);
                            else if(j == 3 + numPartitions && (numTags == 4 + numPartitions))
                                parent = tag;
                            else if(j == 3 + numPartitions &&
                                    (numTags == 5 + numPartitions)) {
                                dom1 = tag;
                                j++;
                                if(fscanf(fp, "%d", &dom2) != 1) {
                                    fclose(fp);
                                    return 0;
                                }
                            }
                        }
                        if(!(numVertices = MElement::getInfoMSH(type))) {
                            if(type != MSH_POLYG_ && type != MSH_POLYH_ &&
                                    type != MSH_POLYG_B) {
                                fclose(fp);
                                return 0;
                            }
                            if(fscanf(fp, "%d", &numVertices) != 1) {
                                fclose(fp);
                                return 0;
                            }
                        }
                    }
                    int *indices = new int[numVertices];
                    for(int j = 0; j < numVertices; j++) {
                        if(fscanf(fp, "%d", &indices[j]) != 1) {
                            delete[] indices;
                            fclose(fp);
                            return 0;
                        }
                    }
                    std::vector<MVertex *> vertices;
                    if(vertexVector.size()) {
                        if(!getMeshVertices(numVertices, indices, vertexVector, vertices,
                                            minVertex)) {
                            delete[] indices;
                            fclose(fp);
                            return 0;
                        }
                    }
                    else {
                        if(!getMeshVertices(numVertices, indices, vertexMap, vertices)) {
                            delete[] indices;
                            fclose(fp);
                            return 0;
                        }
                    }
                    MElement *p = nullptr;
                    bool own = false;

                    // search parent element
                    if(parent != 0) {
                        std::map<int, MElement *>::iterator ite = elems.find(parent);
                        if(ite == elems.end())
                            Msg::Error(
                                        "Parent element (ascii) %d not found for element %d of type %d",
                                        parent, num, type);
                        else {
                            p = ite->second;
                            parents[parent] = p;
                        }
                        std::set<MElement *>::iterator itpo = parentsOwned.find(p);
                        if(itpo == parentsOwned.end()) {
                            own = true;
                            parentsOwned.insert(p);
                        }
                        assert(p != nullptr);
                    }

                    // search domains
                    MElement *doms[2] = {nullptr, nullptr};
                    if(dom1) {
                        std::map<int, MElement *>::iterator ite = elems.find(dom1);
                        if(ite == elems.end())
                            Msg::Error("Domain element %d not found for element %d", dom1,
                                       num);
                        else
                            doms[0] = ite->second;

                        ite = elems.find(dom2);
                        if(ite != elems.end()) doms[1] = ite->second;

                        if(!doms[0])
                            Msg::Error("Domain element %d not found for element %d", dom1,
                                       num);
                        if(dom2 && !doms[1])
                            Msg::Error("Domain element %d not found for element %d", dom2,
                                       num);
                    }
                    delete[] indices;

                    if(elementary < 0) continue;
                    MElement *e = createElementMSH2(this, num, type, physical, elementary,
                                                    partition, vertices, elements,
                                                    physicals, own, p, doms[0], doms[1]);
                    elems[num] = e;
                    elemreg[num] = elementary;
                    elemphy[num] = physical;
                    for(std::size_t j = 0; j < ghosts.size(); j++)
                        _ghostCells.insert(std::pair<MElement *, short>(e, ghosts[j]));
                    if(numElements > 100000)
                        Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
                }
            }
            else {
                int numElementsPartial = 0;
                while(numElementsPartial < numElements) {
                    int header[3];
                    if(fread(header, sizeof(int), 3, fp) != 3) {
                        fclose(fp);
                        return 0;
                    }
                    if(swap) SwapBytes((char *)header, sizeof(int), 3);
                    int type = header[0];
                    int numElms = header[1];
                    int numTags = header[2];
                    int numVertices = MElement::getInfoMSH(type);
                    unsigned int n = 1 + numTags + numVertices;
                    int *data = new int[n];
                    for(int i = 0; i < numElms; i++) {
                        if(fread(data, sizeof(int), n, fp) != n) {
                            delete[] data;
                            fclose(fp);
                            return 0;
                        }
                        if(swap) SwapBytes((char *)data, sizeof(int), n);
                        int num = data[0];
                        int physical = (numTags > 0) ? data[1] : 0;
                        int elementary = (numTags > 1) ? data[2] : 0;
                        int numPartitions = (version >= 2.2 && numTags > 3) ? data[3] : 0;
                        int partition = (version < 2.2 && numTags > 2) ?
                                    data[3] :
                                (version >= 2.2 && numTags > 3) ? data[4] : 0;
                        int parent = (version < 2.2 && numTags > 3) ||
                                (version >= 2.2 && numPartitions &&
                                 numTags > 3 + numPartitions) ||
                                (version >= 2.2 && !numPartitions && numTags > 2) ?
                                    data[numTags] :
                                    0;
                        int *indices = &data[numTags + 1];
                        std::vector<MVertex *> vertices;
                        if(vertexVector.size()) {
                            if(!getMeshVertices(numVertices, indices, vertexVector, vertices,
                                                minVertex)) {
                                delete[] data;
                                fclose(fp);
                                return 0;
                            }
                        }
                        else {
                            if(!getMeshVertices(numVertices, indices, vertexMap, vertices)) {
                                delete[] data;
                                fclose(fp);
                                return 0;
                            }
                        }
                        MElement *p = nullptr;
                        bool own = false;
                        if(parent) {
                            std::map<int, MElement *>::iterator ite = elems.find(parent);
                            if(ite == elems.end())
                                Msg::Error(
                                            "Parent (binary) element %d not found for element %d", parent,
                                            num);
                            else {
                                p = ite->second;
                                parents[parent] = p;
                            }
                            std::set<MElement *>::iterator itpo = parentsOwned.find(p);
                            if(itpo == parentsOwned.end()) {
                                own = true;
                                parentsOwned.insert(p);
                            }
                            assert(p != nullptr);
                        }
                        MElement *e = createElementMSH2(this, num, type, physical,
                                                        elementary, partition, vertices,
                                                        elements, physicals, own, p);
                        elems[num] = e;
                        elemreg[num] = elementary;
                        elemphy[num] = physical;
                        if(numPartitions > 1)
                            for(int j = 0; j < numPartitions - 1; j++)
                                _ghostCells.insert(
                                            std::pair<MElement *, short>(e, -data[5 + j]));
                        if(numElements > 100000)
                            Msg::ProgressMeter(numElementsPartial + i + 1, numElements, true,
                                               "Reading elements");
                    }
                    delete[] data;
                    numElementsPartial += numElms;
                }
            }
            for(int i = 0; i < 10; i++) elements[i].clear();

            std::map<int, MElement *>::iterator ite;
            for(ite = elems.begin(); ite != elems.end(); ite++) {
                int num = ite->first;
                MElement *e = ite->second;
                if(parents.find(num) == parents.end()) {
                    int reg;
                    if(CTX::instance()->mesh.switchElementTags) {
                        reg = elemphy[num];
                    }
                    else {
                        reg = elemreg[num];
                    }
                    switch(e->getType()) {
                    case TYPE_PNT: elements[0][reg].push_back(e); break;
                    case TYPE_LIN: elements[1][reg].push_back(e); break;
                    case TYPE_TRI: elements[2][reg].push_back(e); break;
                    case TYPE_QUA: elements[3][reg].push_back(e); break;
                    case TYPE_TET: elements[4][reg].push_back(e); break;
                    case TYPE_HEX: elements[5][reg].push_back(e); break;
                    case TYPE_PRI: elements[6][reg].push_back(e); break;
                    case TYPE_PYR: elements[7][reg].push_back(e); break;
                    case TYPE_POLYG: elements[8][reg].push_back(e); break;
                    case TYPE_POLYH: elements[9][reg].push_back(e); break;
                    default:
                        Msg::Error("Wrong type of element");
                        fclose(fp);
                        return 0;
                    }
                }
            }
        }
        else if(!strncmp(&str[1], "NodeData", 8)) {
            // there's some nodal post-processing data to read later on, so
            // cache the vertex indexing data
            if(vertexVector.size())
                _vertexVectorCache = vertexVector;
            else
                _vertexMapCache = vertexMap;
            postpro = true;
            break;
        }
        else if(!strncmp(&str[1], "ElementData", 11) ||
                !strncmp(&str[1], "ElementNodeData", 15)) {
            // there's some element post-processing data to read later on
            postpro = true;
            break;
        }

        do {
            if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
        } while(str[0] != '$');
    }

    // store the elements in their associated elementary entity. If the
    // entity does not exist, create a new (discrete) one.
    for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
        _storeElementsInEntities(elements[i]);

    // associate the correct geometrical entity with each mesh vertex
    _associateEntityWithMeshVertices();

    // store the vertices in their associated geometrical entity
    if(vertexVector.size())
        _storeVerticesInEntities(vertexVector);
    else
        _storeVerticesInEntities(vertexMap);

    // store the physical tags
    for(int i = 0; i < 4; i++) _storePhysicalTagsInEntities(i, physicals[i]);

    // copying periodic information from the mesh

    if(!CTX::instance()->mesh.ignorePeriodicity) {
        rewind(fp);

        while(1) {
            while(str[0] != '$') {
                if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
            }
            if(feof(fp)) break;

            if(!strncmp(&str[1], "Periodic", 8) &&
                    strncmp(&str[1], "PeriodicNodes", 13)) {
                readMSHPeriodicNodes(fp, this);
                break;
            }
            do {
                if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
            } while(str[0] != '$');
        }
    }

    fclose(fp);

    if(!CTX::instance()->mesh.ignorePeriodicity) alignPeriodicBoundaries();

    return postpro ? 2 : 1;
}

template <class T>
static void writeElementMSH(FILE *fp, GModel *model, GEntity *ge, T *ele,
                            bool saveAll, double version, bool binary, int &num,
                            int elementary, std::vector<int> &physicals,
                            int parentNum = 0, int dom1Num = 0, int dom2Num = 0)
{
    if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
            && ge->getParentEntity()->dim() > ge->dim())
        return; // ignore partition boundaries

    std::vector<short> ghosts;
    if(model->getGhostCells().size()) {
        std::pair<std::multimap<MElement *, short>::iterator,
                std::multimap<MElement *, short>::iterator>
                itp = model->getGhostCells().equal_range(ele);
        for(std::multimap<MElement *, short>::iterator it = itp.first;
            it != itp.second; it++)
            ghosts.push_back(it->second);
    }

    if(saveAll)
        ele->writeMSH2(fp, version, binary, ++num, elementary, 0, parentNum,
                       dom1Num, dom2Num, &ghosts);
    else {
        if(parentNum) parentNum = parentNum - physicals.size() + 1;
        for(std::size_t j = 0; j < physicals.size(); j++) {
            ele->writeMSH2(fp, version, binary, ++num, elementary, physicals[j],
                           parentNum, dom1Num, dom2Num, &ghosts);
            if(parentNum) parentNum++;
        }
    }

    model->setMeshElementIndex(ele, num); // should really be a multimap...

    if(CTX::instance()->mesh.saveTri && ele->getNumChildren())
        num += ele->getNumChildren() - 1;
}

template <class T>
static void writeElementsMSH(FILE *fp, GModel *model, GEntity *ge,
                             std::vector<T *> &ele,
                             bool saveAll, int saveSinglePartition,
                             double version, bool binary, int &num,
                             int elementary, std::vector<int> &physicals)
{
    if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
            && ge->getParentEntity()->dim() > ge->dim())
        return; // ignore partition boundaries

    // Hack to save each partition as a separate physical entity
    if(saveSinglePartition < 0) {
        if(ele.empty()) return;
        int offset = -saveSinglePartition;
        int maxPhysical = model->getMaxPhysicalNumber(ele[0]->getDim());
        for(std::size_t i = 0; i < ele.size(); i++) {
            std::vector<int> newPhysicals;
            int newElementary = elementary;
            if(elementary > 0) { // classical entity
                newElementary = elementary * offset + ele[i]->getPartition();
                for(std::size_t j = 0; j < physicals.size(); j++)
                    newPhysicals.push_back(physicals[j] * offset +
                                           ele[i]->getPartition());
            }
            else if(elementary < 0) { // partition boundary
                newPhysicals.push_back((maxPhysical - elementary) * offset);
            }
            ele[i]->setPartition(0);
            writeElementMSH(fp, model, ge, ele[i], saveAll, version, binary, num,
                            newElementary, newPhysicals);
        }
        return;
    }

    for(std::size_t i = 0; i < ele.size(); i++) {
        if(saveSinglePartition && ele[i]->getPartition() != saveSinglePartition)
            continue;
        if(ele[i]->getDomain(0)) continue;
        int parentNum = 0;
        MElement *parent = ele[i]->getParent();
        if(parent) parentNum = model->getMeshElementIndex(parent);
        writeElementMSH(fp, model, ge, ele[i], saveAll, version, binary, num,
                        elementary, physicals, parentNum);
    }
}

static int getNumElementsMSH(GEntity *ge, bool saveAll, int saveSinglePartition)
{
    if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
            && ge->getParentEntity()->dim() > ge->dim())
        return 0; // ignore partition boundaries

    int n = 0, p = saveAll ? 1 : ge->physicals.size();

    if(saveSinglePartition < 0 && ge->tag() < 0) p = 1; // partition boundary

    if(saveSinglePartition <= 0)
        n = p * ge->getNumMeshElements();
    else
        for(std::size_t i = 0; i < ge->getNumMeshElements(); i++)
            if(ge->getMeshElement(i)->getPartition() == saveSinglePartition) n += p;
    return n;
}

static int getNumElementsMSH(GModel *m, bool saveAll, int saveSinglePartition)
{
    int n = 0;
    for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(!CTX::instance()->mesh.saveTri) {
            for(std::size_t i = 0; i < (*it)->points.size(); i++)
                if((*it)->points[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
        }
    }
    for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(!CTX::instance()->mesh.saveTri) {
            for(std::size_t i = 0; i < (*it)->lines.size(); i++)
                if((*it)->lines[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
        }
    }
    for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(CTX::instance()->mesh.saveTri) {
            for(std::size_t i = 0; i < (*it)->polygons.size(); i++) {
                int nbC = (*it)->polygons[i]->getNumChildren() - 1;
                n += (saveAll ? nbC : nbC * (*it)->physicals.size());
            }
        }
        else {
            for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
                if((*it)->triangles[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
            for(std::size_t i = 0; i < (*it)->polygons.size(); i++)
                if((*it)->polygons[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
        }
    }
    for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(CTX::instance()->mesh.saveTri) {
            for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++) {
                int nbC = (*it)->polyhedra[i]->getNumChildren() - 1;
                n += (saveAll ? nbC : nbC * (*it)->physicals.size());
            }
        }
        else {
            for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++)
                if((*it)->tetrahedra[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
            for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++)
                if((*it)->polyhedra[i]->ownsParent())
                    n += (saveAll ? 1 : (*it)->physicals.size());
        }
        n -= (*it)->trihedra.size();
    }
    return n;
}

static int _getElementary(GEntity *ge)
{
    if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
            && ge->getParentEntity()->dim() == ge->dim()) {
        // hack for backward compatibility of partitioned meshes in MSH2 format: use
        // elementary tag of parent entity if they are of the same dimension
        // (i.e. if they are not partition boundaries)
        return ge->getParentEntity()->tag();
    }
    else {
        return ge->tag();
    }
}

int GModel::_writeMSH2(const std::string &name, double version, bool binary,
                       bool saveAll, bool saveParametric, double scalingFactor,
                       int elementStartNum, int saveSinglePartition,
                       bool append, bool renumberVertices)
{
    FILE *fp;
    if(append)
        fp = Fopen(name.c_str(), binary ? "ab" : "a");
    else
        fp = Fopen(name.c_str(), binary ? "wb" : "w");
    if(!fp) {
        Msg::Error("Unable to open file '%s'", name.c_str());
        return 0;
    }

    // binary format exists only in version 2
    if(version > 1 || binary)
        version = 2.2;
    else
        version = 1.0;

    // if there are no physicals we save all the elements
    if(noPhysicalGroups()) saveAll = true;

    // get the number of vertices and index the vertices in a continuous
    // sequence

    int numVertices =
            indexMeshVertices(saveAll, saveSinglePartition, renumberVertices);

    // get the number of elements we need to save
    int numElements = getNumElementsMSH(this, saveAll, saveSinglePartition);

    if(version >= 2.0) {
        fprintf(fp, "$MeshFormat\n");
        fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
        if(binary) {
            int one = 1;
            fwrite(&one, sizeof(int), 1, fp);
            fprintf(fp, "\n");
        }
        fprintf(fp, "$EndMeshFormat\n");

        if(numPhysicalNames()) {
            fprintf(fp, "$PhysicalNames\n");
            fprintf(fp, "%d\n", numPhysicalNames());
            for(piter it = firstPhysicalName(); it != lastPhysicalName(); it++) {
                std::string name = it->second;
                if(name.size() > 128) name.resize(128);
                fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
                        name.c_str());
            }
            fprintf(fp, "$EndPhysicalNames\n");
        }

        if(CTX::instance()->mesh.saveTopology) writeMSHEntities(fp, this);

        if(saveParametric)
            fprintf(fp, "$ParametricNodes\n");
        else
            fprintf(fp, "$Nodes\n");
    }
    else
        fprintf(fp, "$NOD\n");

    fprintf(fp, "%d\n", numVertices);

    std::vector<GEntity *> entities;
    getEntities(entities);
    for(std::size_t i = 0; i < entities.size(); i++)
        for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
            entities[i]->mesh_vertices[j]->writeMSH2(fp, binary, saveParametric,
                                                     scalingFactor);

    if(binary) fprintf(fp, "\n");

    if(version >= 2.0) {
        if(saveParametric)
            fprintf(fp, "$EndParametricNodes\n");
        else
            fprintf(fp, "$EndNodes\n");
        fprintf(fp, "$Elements\n");
    }
    else {
        fprintf(fp, "$ENDNOD\n");
        fprintf(fp, "$ELM\n");
    }

    fprintf(fp, "%d\n", numElements);
    int num = elementStartNum;

    _elementIndexCache.clear();

    // parents
    if(!CTX::instance()->mesh.saveTri) {
        for(viter it = firstVertex(); it != lastVertex(); ++it) {
            for(std::size_t i = 0; i < (*it)->points.size(); i++)
                if((*it)->points[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->points[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
        for(eiter it = firstEdge(); it != lastEdge(); ++it) {
            for(std::size_t i = 0; i < (*it)->lines.size(); i++)
                if((*it)->lines[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->lines[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
        for(fiter it = firstFace(); it != lastFace(); ++it) {
            for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
                if((*it)->triangles[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->triangles[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
        for(riter it = firstRegion(); it != lastRegion(); ++it) {
            for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++)
                if((*it)->tetrahedra[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->tetrahedra[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
        for(fiter it = firstFace(); it != lastFace(); ++it) {
            for(std::size_t i = 0; i < (*it)->polygons.size(); i++)
                if((*it)->polygons[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->polygons[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
        for(riter it = firstRegion(); it != lastRegion(); ++it) {
            for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++)
                if((*it)->polyhedra[i]->ownsParent())
                    writeElementMSH(fp, this, *it, (*it)->polyhedra[i]->getParent(), saveAll,
                                    version, binary, num, _getElementary(*it),
                                    (*it)->physicals);
        }
    }
    // points
    for(viter it = firstVertex(); it != lastVertex(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->points, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // lines
    for(eiter it = firstEdge(); it != lastEdge(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->lines, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // triangles
    for(fiter it = firstFace(); it != lastFace(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->triangles, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // quads
    for(fiter it = firstFace(); it != lastFace(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->quadrangles, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // polygons
    for(fiter it = firstFace(); it != lastFace(); it++) {
        writeElementsMSH(fp, this, *it, (*it)->polygons, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // tets
    for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->tetrahedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // hexas
    for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->hexahedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // prisms
    for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->prisms, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // pyramids
    for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->pyramids, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // polyhedra
    for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->polyhedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
    }
    // level set faces
    for(fiter it = firstFace(); it != lastFace(); ++it) {
        for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
            MTriangle *t = (*it)->triangles[i];
            if(t->getDomain(0))
                writeElementMSH(fp, this, *it, t, saveAll, version, binary, num,
                                _getElementary(*it), (*it)->physicals, 0,
                                getMeshElementIndex(t->getDomain(0)),
                                getMeshElementIndex(t->getDomain(1)));
        }
        for(std::size_t i = 0; i < (*it)->polygons.size(); i++) {
            MPolygon *p = (*it)->polygons[i];
            if(p->getDomain(0))
                writeElementMSH(fp, this, *it, p, saveAll, version, binary, num,
                                _getElementary(*it), (*it)->physicals, 0,
                                getMeshElementIndex(p->getDomain(0)),
                                getMeshElementIndex(p->getDomain(1)));
        }
    }
    // level set lines
    for(eiter it = firstEdge(); it != lastEdge(); ++it) {
        for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
            MLine *l = (*it)->lines[i];
            if(l->getDomain(0))
                writeElementMSH(fp, this, *it, l, saveAll, version, binary, num,
                                _getElementary(*it), (*it)->physicals, 0,
                                getMeshElementIndex(l->getDomain(0)),
                                getMeshElementIndex(l->getDomain(1)));
        }
    }

    if(binary) fprintf(fp, "\n");

    if(version >= 2.0) {
        fprintf(fp, "$EndElements\n");
    }
    else {
        fprintf(fp, "$ENDELM\n");
    }

    writeMSHPeriodicNodes(fp, entities, renumberVertices);

    fclose(fp);

    return 1;
}

int GModel::_writePartitionedMSH2(const std::string &baseName, bool binary,
                                  bool saveAll, bool saveParametric,
                                  double scalingFactor)
{
    int numElements;
    int startNum = 0;
    for(std::size_t partition = 1; partition <= getNumPartitions();
        partition++) {
        std::ostringstream sstream;
        sstream << baseName << "_" << partition << ".msh";

        numElements = getNumElementsMSH(this, saveAll, partition);
        Msg::Info("Writing partition %d in file '%s'", partition,
                  sstream.str().c_str());
        _writeMSH2(sstream.str(), 2.2, binary, saveAll, saveParametric,
                   scalingFactor, startNum, partition, false,
                   false); // NO RENUMBERING!
        startNum += numElements; // update for next iteration in the loop
    }

#if 0
    if(_ghostCells.size()){
        Msg::Info("Writing ghost cells in debug file 'ghosts.pos'");
        FILE *fp = Fopen("ghosts.pos", "w");
        fprintf(fp, "View \"ghosts\"{\n");
        for(std::multimap<MElement*, short>::iterator it = _ghostCells.begin();
            it != _ghostCells.end(); it++)
            it->first->writePOS(fp, false, true, false, false, false, false);
        fprintf(fp, "};\n");
        fclose(fp);
    }
#endif

    return 1;
}
